Prom antiferromagnetism to d-wave superconductivity in the 2D t — J model 
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We have found that the two dimensional t-J model, for the physical parameter range J/t = 0.4 
reproduces the main experimental qualitative features of High-Tc copper oxide superconductors: 
d-wave superconducting correlations are strongly enhanced upon small doping and clear evidence 
of off diagonal long range order is found at the optimal doping S ~ 0.15. On the other hand 
antiferromagnetic long range order, clearly present at zero hole doping, is suppressed at small hole 
density with clear absence of antiferromagnetism at 5 >'~ 0.1. 
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The interplay between antiferromagnetism and super- 
conductivity in the Cu02 layers of the high Tc com- 
pounds is one of the most important effects where strong 
electron correlation may play the main role. How- 
ever, after many years of theoretical studies and exper- 
imental efforts [H the most obvious question is still un- 
clear: whether the occurrence of high Tc superconduc- 
tivity is determined by the proximity of the compound 
to a perfect antiferromagnetic insulator. 

In case strong correlation is the dominant force driv- 
ing from antiferromagnetism to superconductivity a well 
accepted model is the 2D t-J model [|j: 



H ^ J {S^■ Sj- ^n^Uj) - t i^^Cja + h.c), 



<i,J> 



<h3> 



(1) 



where c]^ creates an electron of spin a on the site i, Ui and 
Si being the electron number and spin operators respec- 
tively. Double occupations are forbidden and < i,j > 
denotes nearest neighbor summation over the L lattice 
sites with periodic boundary conditions. 

In the last decade the investigation of the properties of 
the 2D t — J model (and of the parent Hubbard model) 
has been a challenge for numerical calculations. Exact di- 
agonalization (ED) [Q shows that antiferromagnetic cor- 
relations are resistant up to (5 ~ 0.15 and superconductiv- 
ity is present at intermediate doping but the lattice sizes 
considered were too small for being conclusive. On the 
contrary the quantum Monte Carlo (QMC) methods al- 
low simulations on larger systems but suffer from the well 
known "minus sign problem" instability, which makes the 
simulation impossible at low enough temperatures. 

At present, this instability can be controlled, only at 
the price of introducing some approximation, such as 
the fixed node (FN) approximation which is strictly 
variational on the ground state energy, the constrained 
path quantum Monte Carlo ^ and the Green function 
Monte Carlo with stochastic reconfiguration GFMCSR 
0, which has been developed to improve the accuracy 
of the FN. Both the FN and GFMCSR techniques will 
be extensively used in this work. Similar approximations 



on the ground state wavefunction can be obtained by 
applying one (or more) Lanczos steps (LS) to the varia- 
tional wavefunction , or also using the density ma- 
trix renormalization group (DMRG), which in 2D is also 
affected by a sizable error, and is not "numerically exact" 
as in ID. 

All these approximations allow to obtain a rather accu- 
rate value of the ground state energy of the model, with 
an error typically less than 1% of the correlation energy 
even for large L. However this kind of accuracy for the 
energy certainly does not allow to draw reasonable con- 
clusions on the interesting long range properties of the 
model, see e.g. [|ll|. On the other hand it is reasonable 
to expect that, by using approximate techniques that do 
not spoil the local character of the Hamiltonian, a similar 
good accuracy can be obtained on the ground state ex- 
pectation value of short range operators like, for instance, 
the kinetic energy and the exchange energies in Eqj^. 
These operators O, acting only on nearest neighbor sites, 
share the important property that, if added to the Hamil- 
tonian ( Hh H — hO ) do not change its local charac- 
ter. Moreover this kind of perturbation typically leads to 
a sizable change of the ground state energy per site Eh 
even in the linear regime Eh = Eq — h < O > /L + o{h) , 
providing a very reliable estimate of the ground state 
expectation value < O >, as the energy E{h) can be 
accurately determined for few values of the field h. 

So far, in the literature |12jl3,0], the ground state 



expectation value of the squared order parameter O is 
estimated on an approximate ground state ipo, by taking, 
simply, its bare expectation value < ^o\0'^\tpo >. For 
long range operators such as O^, this may lead to very 
poor approximations, unless the method is almost exact. 

In order to detect superconducting long range order 
with a more controlled approximation, we perform simu- 
lations in the grand canonical ensemble and add to iJ a 
short range operator which creates or destroys a d-wave 
singlet Cooper pair: 



H{h) = H -h{A- 



(2) 



where A+ = E<^,J> M,jic+c+^ + c+c+) and My = 1 
or —1 if the bond < i,j > is in the x or y direction re- 
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spectively, while /z is the chemical potential and N the 
particle operator. FN and GFMCSR allow to compute 
quite accurately the ground state energy E{h) also in 
presence of the field h. To this purpose a fundamental 
role is played by the guiding wavefunction which allows to 
perform importance sampling. We generalize the A'^ par- 
ticle, d-wave symmetry, BCS guiding wavefunction 
{\BCS)) to the grand canonical ensemble by introducing 
a proper weight /jv for each N particle sector: 



(3) 



N 



where Va projects out doubly occupied sites and Pjv se- 
lects the A^— particle component of the wavefunction. 
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FIG. 1. Pi for = 16, L = 18 

Our purpose is to compute the anomalous average of 
the order parameter = K^V + 2|A+|7V)|/L, where \N) 
and -I- 2) are the N and -I- 2 particle ground state 
respectively, pd can be non zero even on a finite size and 
zero external field. Moreover, if superconducting long 
range order occurs, remains finite for L — > cxd. 

In order to compute on finite size systems we have 
implemented the following simple strategy. We choose 
the chemical potential /i in a way that the ground state 
energies per site Ey; and £'Ar+2 for the A^ and A^ -I- 2 
particles are degenerate. In order to reduce the ground 
state energy statistical error we optimize the variational 
parameters /jv by restricting ourselves to the subspaces 
of A^ and A^ -I- 2 particles relevant for the matrix element 
Vd, !n being zero otherwise. In the guiding function 
/tv and Jn+i are then determined by requiring that the 
average particle number < '\1)q\N\iI}q > is equal to A^ -I- 
1. The first order correction to the energy due to the 
perturbation (j^) in this restricted Fock space is given by 
the eigenvalues of the secular matrix: 



En ±hpd 
±hpd En+2- 



(4) 



It easily follows that E{h) ~ Em — hpd, meaning that 
the anomalous average of the order parameter can be 
computed as an energy difference (E^ — E(h))/h for h 
0. A long range property of the model can be probed by 



studying the ground state energy change under the effect 
of a local perturbation. We expect this scheme to be a 
much more controlled and accurate way to characterize 
the long distance behavior of a model. 
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FIG. 2. VMC(arrows), FN(empty dots) and GFMCSR(full 
dots) calculation of at J = 0.4t. (see text for details) 

As can be seen from the comparison with the exact 
results in Fig. 0(a), at J = 0.4i, the VMC highly overes- 
timates the order parameter. The FN reduces this value. 
The GFMCSR, implemented by reconfiguring the unper- 
turbed energy of the two subspaces at A^ and A^ -I- 2 par- 
ticles in an independent way, is almost exact. 

In order to attempt a finite size scaling for the order 
parameter we compute pd for much larger sizes (Fig. |^). 
As can be seen in the i = 50 lattice case both the FN 
and the GFMCSR reduce the variational value. We have 
tested the accuracy of the calculation and the depen- 
dency of the results from the chosen guiding wavefunc- 
tion, by reducing the optimal energy variational parame- 
ter Adw = 0.65 (dots connected by full lines in Fig. ^ a) 
to the value of ^dw — 0.3 (dashed lines). This implies a 
sizable reduction of pd within VMC. The FN evaluation 
of Pd correctly enhances this value, getting closer to the 
more reliable estimate obtained with the optimal energy 
variational parameter. The GFMCSR method, the most 
accurate technique used here, is, remarkably, rather in- 
sensitive to the choice of the guiding function, being the 
difference for the two GFMCSR results a conservative es- 
timate of the possible error in the determination of pd- 
GFMCSR seems to improve by the same amount the FN 
estimate of pd, both for the 18 sites (fig. |l|) and 50 sites 
(fig. ^ a), and this improvement is expected to remain 
even for larger sizes, being GFMCSR, as well as FN, a 
size consistent approximation. 

The 98-site calculation shows that the VMC value of 
Pd is enhanced both by the FN and GFMCSR calculation 
and remarkably the computed value is very close to the 
one obtained for the 50 site lattice. 

Our results at this doping and J jt value display , all 
consistently, stronger and stronger d-wave correlations, 
as the accuracy of our numerical techniques are improved 
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and lattice size increased. We believe that this represents 
a robust evidence of d-wave superconductivity in the 2D 
t — J model. However the limited number of lattice sizes 
considered does not allow us to perform an accurate finite 
size scaling. As shown in Fig. ^, size effects are present 
also at the variational level and the true order parameter 
maybe below the value ^ 0.12 reported in the picture. 
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FIG. 3. Size scaling of pd sX J = 0.4t. Lines connecting FN 
and GFMCSR in (a) are guides to the eye, least square fit for 
the variational method in (b). 

Since the t — J model originates from the doping of 
an antiferromagnetic Mott insulator it is interesting to 
understand if the antiferromagnetic character of the un- 
doped ground state is resistant upon doping. Follow- 
ing a similar procedure to the one used for the super- 
conducting long range order we added to the Hamil- 
tonian a short range perturbation coupled to the stag- 
gered magnetizationim^ = l/^X^i? ; namely 
H ^ H — IT'^rSr{~)^i ^iid compute rrih in presence 
of the field h either by differentiating the energy per site 
ruh = —dE(h)/dh or by using the forward walking tech- 
nique, whenever possible (FN) . ||l5| For this quantity the 
FN and GFMCSR are consistent for small field, meaning 
that the FN is already enough accurate for the magnetic 
phase diagram. 

For the Heisenberg antiferromagnet, where broken 
symmetry occurs, the magnetization as a function of the 
rescaled field /i — > /i • L lies on a universal curve fl^ j 
which weakly depends on the system size. This size de- 
pendence is almost negligible if compared to the one af- 
fecting the squared order parameter (Fig. || b) ||l5|,|l^ . 
This feature strengthen the validity of our results that are 
all based upon ground state expectation values of short 
range operators in presence of a field. At finite doping, 
computationally heavier, we have chosen to work with 
a single field for each size and tuned at zero doping in 
order to reproduce on the available finite systems the in- 
finite size order parameter: h = x/ JL with x = 0.392. 
It turns out that the antiferromagnetic correlations are 
present even at finite doping up to 5c — 0.10 see Fig.H 



(b), in qualitative agreement with experimental findings 
(J^^p ^ 3 — 5%). For a quantitative agreement, other 
terms must be probably added to the Hamiltonian, as 
suggested in |jl^ . In the optimal doping region the stag- 
gered magnetization is vanishingly small even in presence 
of a sizable magnetic field, meaning that long range order 
has disappeared favoring a pure d — wave superconduct- 
ing state. 
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FIG. 4. Staggered magnetization for x = hJL — 0.392 
(a), nih for 5 — (b). The horizontal dotted line repre- 
sents the squared order parameter value. Remaining lines are 
guides to the eye. 

The interplay between antiferromagnetism and super- 
conductivity appears to be a fundamental point in the 
phase diagram of the t — J model. For small doping the 
matrix element < A^-|-2|A"'"|A'' >, is strongly suppressed 
but antiferromagnetism is still present. Close to the Mott 
insulator ((5 = 0), as pointed out previously |ll9|-^, there 
is a strong tendency to have a phase separation instabil- 
ity between an hole rich uniform phase and an undoped 
antiferromagnetic insulator. In the phase diagram shown 
in Fig. (|^) , that we have obtained with the same method 
(GFMCSR using only the energy in the reconfiguration 
scheme) used for the computation of the d-wave order pa- 
rameter, the PS boundary is quite far from the optimal 
doping region at J/t — 0.4. However the compressibility 
of the electron system is very large (^ w 0.54t) about 
20 times larger than the corresponding spinless fermion 
compressibility, in surprising numerical agreement with 
a spinless fermion model with renormalized flat band. 
p3[ Thus the proximity to an antiferromagnetic insu- 
lator strongly enhances charge fluctuations determining 
-for physical J/ 1 values- a d-wave superconducting phase 
before the PS instability. 

We believe that a large value of the compressibility 
is very important to stabilize superconductivity even 
in presence of long range Coulomb repulsion, certainly 
present in the physical system but missing in the t — J 
model. For large compressibility the Thomas-Fermi 
screening length f = 2ire^Hn lUl very short compared 
to the Cu-Cu distance, so that the screening is very much 
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effective. We have verified this picture (on smaller sizes) 
by adding to the t — J model a repulsive nearest-neighbor 
interaction ^X)<ij>"'»% found still strong super- 
conducting correlations, weakly suppressed even for large 
VI J - 1. 

Another mechanism in competition with superconduc- 
tivity, is the formation of so called "stripes" in the ground 
state of the t — J model, as recently found by White and 
Scalapino with DMRG. |l3| In order to test this pos- 
sibility we have compared our results with the DMRG 
ones on a 12 X 6 system with 8 holes and open boundary 
conditions at J/t = 0.4. In this case, DMRG is quite 
more accurate than our techniques in the energy esti- 
mate, but it is not yet clear whether the same is true 
for correlation functions especially the ones described by 
short range operators, where our approximations seem 
to be quite reliable. Within this accuracy for correla- 
tion functions, we have not found any clear indication of 
"stripes", in qualitative disagreement with the DMRG 
results, and confirming our previous work [2^ ], obtained 
with periodic boundary conditions. In this case, remark- 
ably, the possibility to use translation invariance, allows 
calculations by far more accurate and reliable even com- 
pared with the best DMRG results, both for energies and 
correlation functions. 



sizable systematic errors on correlation functions. With 
our technique we overcome this difficulty by estimating 
only short range operators expectation values with en- 
ergy difference calculations. The short range operators 
expectation values are less sensitive to finite size effects, 
and contain the useful information to establish absence 
or presence of long range order. 

We finally remark that it is extremely important to use 
a very accurate method to rule out superconductivity at 
small doping for a strongly correlated system like the t-J 
model. Even at the variational level the superconduct- 
ing order parameter that is very large before Gutzwiller 
projection, becomes an extremely small quantity after 
this projection (see Fig ^ b). This strong suppression of 
d — wave pairing, can be easily shown at the variational 
level (see Fig.||) and proven aX 6 — Q, and is a crucial 
property of strongly correlated systems. We acknowledge 
S. White for sending us numerical results before publi- 
cation, M. Fabrizio, L. Capriotti, M. Capone, F. Becca 
for useful discussion and A. Parola for careful reading of 
the manuscript This work was supported by INFM (PRA 
HTCS), and partly by MURST COFIN97 and COFIN99. 
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FIG. 5. Instability (PS) of the uniform phase evaluated by 
GFMCSR using the Maxwell construction for the 98 site lat- 
tice. Errors are estimates of finite size effects and correspond 
to twice the difference between the 98 and 50 site critical dop- 
ing |2^. SC label the 5 = 0.14 were Pd has been computed, 
AF label the antiferromagnetic region 

In this model we thus recover the most simple scenario 
, appeared in the early days of HTc superconductivity, 
namely that the strong correlation alone may drive the 
system from antiferromagnetism to superconductivity. 

The contradictory results present in the literature so 
far are, in our opinion, mainly due to the general attempt 
of computing a long range quantity by using approxima- 
tions that weakly affect energy estimates but may lead to 
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